Introduction

Depression is one of the most common mental health conditions worldwide. It can have significant consequences for individuals, communities, and wider society, including inmpacts on wellbeing and employment. One form of depression known to worsen during the winter months is Seasonal Affective Disorder (Rosenthal,M.B(2004)).

The Covid-19 pandemic created major distruptions to daily life, including prolonged periods spent indoors due to lockdowns. Reduced exposure to natural light is a known risk factor for SAD(Rosenthal,M.B(2004)), and the combination of winter daylight shortages and enforced indoor isolation may have contributed to increased levels of depressive symptoms during the pandemic. These changes may be reflected in pattern of antidepressant prescribing, particularly during winter periods.

The aim of this study is to investigate wether the COVID-19 pandemic influenced antidepressant prescribing and whether there are inequalities or geographic differences associated with this trend.

Methods

All data used in this study was obtained from publicly available NHS Scotland and Scottish Government sources.Antidepressant prescribing data were downloaded from the Data by Prescriber Location dataset, focusing specifically on winter periods (December to Februray) for the years 2019-2023, from the Public Health Scotland OpenData portal. Winter population estimates were taken from the National Records of Scotland census release.Socioeconomic information was sourced from the Scottish Index of Multiple Deprivation (SIMD 2020v2) datasets published by the Scottish Government. Spatial analysis was supported by incorporating NHS Scotland Health Board boundary shapefiles, which are publicly available through PHS and used to map regional variation in prescribing.

Antidepressant medications were identified using the BNF item code prefix “0403”, and all prescriptions with codes starting with “0403” were retained for analysis. These filtered data were then processed and summarised to generate visualisations of prescribing patterns.

For the purpose of this study, the COVID-19 period is defined as March 2020 to June 2021, corresponding to the main period of national lockdowns. This places Winter 2019 as pre-COVID, Winter 2020/2021 as during COVID, and Winter 2022 as post-COVID.

Figure1,Figure2,Figure3 are interactive - hover over the elements to see additional details.

##Loading Necessary Packages
library(tidyverse) 
library(here) # directory stucture
library(gt) # tables
library(janitor) # cleaning data
library(ggplot2) # plotting graph
library(sf) # to read in map data 
library(readxl) # to read in map data
library(plotly) # to make interactive
library(viridis) # colour palette for map
library(ggiraph) # make ggplot2 graphs interactive
# Loading multiple data files (18 CSVs) containing counts of medications dispensed in the community
files <- list.files(here("data", "winter_data"), pattern = "csv")
winter_data <- files %>% # Read all downloaded winter prescribing files stored in the 'winter_data' folder
  map_dfr(~read_csv(here("data", "winter_data", .))) %>% 
clean_names()#combine into one dataframe and clean column names
filtered_winter_data <- winter_data %>% 
filter(str_starts(bnf_item_code,"0403")) %>%  #Filter for antidepressants only (BNF item code starts with "0403" 
  mutate(year = as.numeric(substr(paid_date_month,1,4)), month = as.numeric(substr(paid_date_month,5,6))) %>% #Extract year and month from 'paid_date_month' for grouping into winter periods
  mutate(winter_year=case_when(month == 12 ~ year + 1, 
month %in% c(1,2) ~ year) )  #Create 'winter_year' column: December counts belong to the next year, Jan-Feb to current year
filtered_winter_data <- filtered_winter_data %>% 
  unite("healthboards",hbt2014,hbt,sep = "_") #Merge Health Board columns into a single column since they contain the same information
filtered_winter_data$healthboards <- gsub("[NA]","",filtered_winter_data$healthboards) 
    filtered_winter_data$healthboards <-
      gsub("_","",filtered_winter_data$healthboards)#Remove any remaining 'NA' strings and underscores to standardize Health Board codes

Figure 1

#summarise total antidepressant prescriptions for each winter year
winter_years_data <- filtered_winter_data %>% 
  group_by(winter_year) %>% 
  summarise(total_items=sum(number_of_paid_items,na.rm = TRUE))
#plot line graph to visualise overall trend in total prescriptions across winter years
plot <- suppressWarnings( ggplot(winter_years_data,aes(x=winter_year,y=total_items)) + #line connecting total prescriptions per year
  geom_line(linewidth=0.7,colour = "blue") +
  geom_point(aes(text = paste0(
      "<b>Winter Year:</b> ",winter_year, "<br>",
      "<b>Total Items:</b> ",scales::comma(total_items)),
size = 1))+ #points with interactive tooltip showing year and total
  geom_vline(xintercept = 2020, 
linetype = "solid", colour = "brown") + geom_vline(xintercept = 2022, 
linetype = "solid", colour = "brown") + #reference lines to mark COVID-related winters
  scale_x_continuous(breaks=2017:2023) + #ensures all years appear on x-axis
  labs(title="Total Antidepressant Prescriptions During Winter Season",x="Winter Year",y="Total Antidepressant Prescriptions") +
  theme_minimal(base_size = 10))#thiscode suppresses warnings
ggplotly(plot,tooltip="text",height=250) #convert ggplot to interactive plotly object with tooltips

The line graph shows a steady upward trend in winter antidepressant prescribing across all years. Although the COVID-19 period (2020–2022) is highlighted, prescribing patterns do not markedly diverge from pre-pandemic trends, with increases during COVID appearing slightly less steep than earlier years. The y-axis does not start at 0 to improve visibility of year-to-year changes. 

#Prepare Population
population <- readxl::read_excel(here("data","population.xlsx"), skip=10) %>% 
  clean_names() %>% 
  select(x2,all_people) %>% 
  filter(!is.na(all_people)) %>% 
  rename(h_bname = "x2",hb_population = "all_people") %>% 
filter(!str_detect(hb_population,"Cells"))

filtered3_winter_data <- filtered_winter_data %>% 
  group_by(healthboards,winter_year,gp_practice) %>% 
  summarise(paid_quantity = sum(paid_quantity,na.rm=TRUE)) 

SIMD <- readxl::read_excel(here("data","SIMD.xlsx")) %>% 
  clean_names() 

combined_pop_simd <- SIMD %>% 
  full_join(filtered3_winter_data,join_by(h_bcode == healthboards)) %>% 
  left_join(population,join_by(h_bname == h_bname))

Figure 2

# Make sure 'combined' exists as in previous chunks

hb_prescribing <- filtered_winter_data %>%
filter(winter_year >= 2019) %>%
group_by(healthboards, winter_year) %>%
summarise(total_paid = sum(paid_quantity, na.rm = TRUE))

simd_clean <- SIMD %>% select(h_bcode, h_bname) %>% distinct()

combined <- hb_prescribing %>%
left_join(simd_clean, join_by(healthboards == h_bcode)) %>%
left_join(population, by = "h_bname") %>%
mutate(quantity_per_head = total_paid / hb_population) %>%
filter(!is.na(h_bname))

# Load shapefile for NHS Health Boards
NHS_healthboards <- st_read(here("data","Week6_NHS_HealthBoards_2019.shp"), quiet = TRUE) %>%
clean_names() %>% rename(h_bname = hb_name)

retry_mapped_data <- NHS_healthboards %>%
left_join(combined, by = "h_bname") %>%
st_as_sf() %>%
mutate(quantity_per_head = replace_na(quantity_per_head, 0))

# Plot interactive map
plot_map <- retry_mapped_data %>%
ggplot() +
geom_sf_interactive(
aes(fill = quantity_per_head,
tooltip = paste0(h_bname,
"\nWinter Year: ", winter_year,
"\nItems per Head: ", round(quantity_per_head,2))),
colour = "white",
size = 0.1) +
scale_fill_viridis_c(option = "mako", direction = -1, name = "Items per Head") +
labs(title="Antidepressant Prescriptions per Head",
subtitle="By Health Board and Winter Year") +
facet_wrap(~winter_year,nrow = 2) +
theme_void(base_size = 8) +
theme(strip.text = element_text(size=10, face="bold"),
plot.title = element_text(face="bold", size=12),
plot.subtitle = element_text(size=8))
# Convert to interactive map
interactive_map <- girafe(ggobj = plot_map, width_svg = 10, height_svg = 6)
interactive_map

The map shows the geographic distribution of antidepressant prescriptions per head, revealing a consistent north–south gradient with higher rates in central urban regions. Post-2020 intensification suggests a system-wide rise in prescribing, with areas such as Lanarkshire and Glasgow showing the largest increases compared with 2019. Regions with historically higher prescribing appear to have experienced the greatest growth during and after the pandemic.

Figure 3

winter_population <- list.files(here("data", "winter_population"), pattern = "csv", full.names = TRUE) %>%
  map_dfr(read_csv) %>%
  clean_names() %>%
  select(date, practice_code, hb, all_ages) %>%
  mutate(winter_year = as.numeric(substr(date, 1, 4)),
    winter_year = case_when(
      winter_year == 2020 ~ 2019,
      winter_year == 2023 ~ 2022,
      TRUE ~ winter_year ),
    healthboards = hb) %>%
  select(-hb)
#summarise and filter combined population data for selected years
summary_winter_with_simd <- combined_pop_simd %>%
  filter(!is.na(simd2020v2_quintile), winter_year %in% c("2019","2022")) %>%
  group_by(simd2020v2_quintile, h_bname, gp_practice, winter_year) %>%
  summarise(avg_paid_over_winters = mean(paid_quantity), .groups = "drop") %>%
  rename(practice_code = gp_practice) %>%
  mutate(winter_year = as.numeric(winter_year))
# Join datasets
joined_data <- summary_winter_with_simd %>%
  full_join(winter_population, by = c("winter_year", "practice_code"))
# Calculate prescriptions per 1000
final_boxplot <- joined_data %>%
  group_by(simd2020v2_quintile, h_bname, practice_code, winter_year) %>%
  reframe(
    avg_per_1000 = (avg_paid_over_winters / all_ages) * 1000,
    year_label = factor(winter_year,
                        levels = c(2019, 2022),
                        labels = c("Pre-COVID (2019)", "Post-COVID (2022)"))
  )

#create boxplot of average prescriptions per 1000 by SIMD quintile, faceted by pre/post COVID
box2 <- final_boxplot %>%
  ggplot(aes(x = simd2020v2_quintile, y = avg_per_1000, fill = factor(simd2020v2_quintile))) +
  geom_boxplot() +
  scale_fill_viridis(discrete = TRUE, alpha = 0.6) +
  scale_y_continuous(labels = scales::label_comma()) +
  coord_cartesian(ylim = c(0, 70000))  + #scaled y axis graph to be able to visualise better
  theme_minimal(base_size=10) +
  theme(legend.position = "none") +
  facet_wrap(~year_label) +
  labs(title = "Average Winter Antidepressant Prescriptions Per by SIMD Quintile",
    x = "SIMD Quintile (1 = Most Deprived)",
    y = "Avg prescriptions per 1000")

ggplotly(box2, width=900,height=500)

The boxplot shows no strong socioeconomic gradient in antidepressant prescribing across SIMD quintiles, with median rates remaining relatively flat in both 2019 and 2022. Although overall prescribing increased uniformly after COVID-19, greater variability and more extreme outliers in 2022 suggest widening differences between individual health boards or practices, rather than between deprivation groups.

Figure 4

#create population-weighted deprivation profile per health board
hb_deprivation <- SIMD %>%
  group_by(h_bname) %>%
  summarise(
#weighted mean SIMD quintile (lower value = more deprived)
    avg_quintile = weighted.mean(simd2020v2_quintile, w = population, na.rm = TRUE),
#percentage of population living in most deprived quintiles (Q1 & Q2)
    pct_deprived = sum(population[simd2020v2_quintile <= 2], na.rm = TRUE) / 
                   sum(population, na.rm = TRUE) * 100,
#total population per health board (for reference)
    total_pop = sum(population, na.rm = TRUE)
  ) %>%
  arrange(avg_quintile)  #rank health boards from most to least deprived
#join deprivation ranking to prescribing data (2019 = pre-COVID, 2022 = post-COVID)
final_lollipop_data <- combined %>% 
  filter(winter_year %in% c("2019", "2022")) %>%
  mutate(period = ifelse(winter_year == 2019, "pre", "post")) %>%
  full_join(hb_deprivation, by = "h_bname")
#calculate percentage change in prescribing per health board
change_summary <- final_lollipop_data %>% 
  group_by(h_bname, period) %>% 
  summarise(avg_prescribe = mean(total_paid),     
    avg_quintile = first(avg_quintile),
    pct_deprived = first(pct_deprived),
    total_pop = first(total_pop)) %>% 
  pivot_wider(
    names_from = period, 
    values_from = avg_prescribe
  ) %>%
  mutate(
    pct_change = round(((post - pre) / pre) * 100, 1),
    avg_quintile = round(avg_quintile, 2),
    pct_deprived = round(pct_deprived, 1),
    pre = round(pre, 0),
    post = round(post, 0)
  ) %>%
  arrange(desc(pct_change))  #sort by percentage change (highest first)

# Create formatted table
prescribing_table <- change_summary %>%
  select(h_bname, pre, post, pct_change, avg_quintile, pct_deprived) %>%
  rename(
    "Health Board" = h_bname,
    "Pre-COVID (2019)" = pre,
    "Post-COVID (2022)" = post,
    "% Change" = pct_change,
    "Avg SIMD Quintile" = avg_quintile,
    "% in Deprived Areas (Q1-Q2)" = pct_deprived
  )

prescribing_table %>% 
gt() %>% 
tab_header(title=md("**Antidepressant Prescribing** Changes by Health Board **Before and After COVID-19** (2019-2022)"),subtitle=md("with *deprivation indicators*")) %>% 
tab_spanner(label="average total prescriptions",columns = c('Pre-COVID (2019)','Post-COVID (2022)')) %>% 
tab_spanner(label="Deprivation",columns = c('Avg SIMD Quintile','% in Deprived Areas (Q1-Q2)')) %>% 
  tab_options(
    table.font.size = px(12),
    data_row.padding = px(2),
    heading.padding = px(2) ) %>% 
   tab_footnote(footnote = md("**Avg SIMD Quintle** 1 represents most deprived areas in scotland ")) %>%
tab_caption(md("The table shows that Glasgow and Lanarkshire had the largest increases in antidepressant prescribing after COVID-19. This cannot be explained by average deprivation alone, but becomes clearer when considering that both regions have the highest concentrations of people living in the most deprived areas. Rural and island boards showed far smaller increases, suggesting that urbanisation and population density played a bigger role than deprivation in driving prescribing changes during COVID."))
The table shows that Glasgow and Lanarkshire had the largest increases in antidepressant prescribing after COVID-19. This cannot be explained by average deprivation alone, but becomes clearer when considering that both regions have the highest concentrations of people living in the most deprived areas. Rural and island boards showed far smaller increases, suggesting that urbanisation and population density played a bigger role than deprivation in driving prescribing changes during COVID.
Antidepressant Prescribing Changes by Health Board Before and After COVID-19 (2019-2022)
with deprivation indicators
average total prescriptions
 % Change
Deprivation
Pre-COVID (2019) Post-COVID (2022) Avg SIMD Quintile % in Deprived Areas (Q1-Q2)
Greater Glasgow and Clyde
6555083 23757849 262.4 2.70 52.0
Lanarkshire
3666290 13186662 259.7 2.67 51.8
Fife
6681894 7813849 16.9 3.03 39.8
Grampian
7764826 9058806 16.7 3.62 20.6
Western Isles
366503 426735 16.4 2.85 15.1
Highland
4424027 5124922 15.8 3.12 25.5
Lothian
13123200 15122300 15.2 3.40 32.0
Orkney
319190 367197 15.0 3.52 16.3
Dumfries and Galloway
2563265 2936457 14.6 2.94 33.9
Ayrshire and Arran
6720774 7619110 13.4 2.60 52.5
Forth Valley
5275843 5949958 12.8 3.14 36.1
Shetland
305882 344729 12.7 3.52 6.0
Tayside
6709598 7491600 11.7 3.11 34.8
Borders
1992572 2211388 11.0 3.23 21.7
Avg SIMD Quintle 1 represents most deprived areas in scotland

Conclusion This analysis shows that antidepressant prescribing increased across almost all Scottish health boards after COVID-19, with the largest rises seen in Glasgow and Lanarkshire. However, average SIMD quintile did not explain these differences, as no clear deprivation gradient appeared before or after the pandemic Figure 3.

When population-weighted deprivation was included, a clearer pattern emerged: health boards with the highest proportions of residents in the most deprived areas (Q1–Q2) experienced the greatest increases Figure 4. This suggests that COVID-19 may have intensified pre-existing mental health inequalities in large, urban, socioeconomically deprived regions.

This aligns with Figure 2, which shows higher prescribing per head in central Scotland’s densely populated urban belt compared with more rural northern areas, highlighting the influence of population density and urbanisation.

Overall, although prescribing did rise during COVID-19, the trend shown in Figure 1 indicates that this increase was part of a longer-term upward trajectory rather than a sharp pandemic-specific surge.

Limitations Several limitations reduce the strength of these findings. Firstly, the data reflects prescribing volume, not whether the antidepressants were newly started or continued, so it may not directly truly represent true changes in mental health across deprivation levels.In addition, important confounding variables such as changes in healthcare access,remote consultations, and shifting service pressures durung COVID-19 potentially affecting the accuracy of the results.Finally, the analysis was only done at health board level, so findings cannot be generalised to individual gp practice levels so GP practices, and no conclusions can be made about patient-level experiences or need.

Next Steps Future work could be done to link prescribing levels to practice level SIMD and demographic characteristics to better understand drivers of change. Work should be done analysing antidepressant prescriptions all year round and compare it with the winter-only dataset used here, to determine whether truly seasonal peaks exist and whether winter peaks changed post-COVID. Additionally, applying time-series modelling would confirm whether the post COVID rise represents a temporary shock or a long term behavioural shift in prescribing patterns.

References National Records of Scotland (2022). Scottish Census. Scottish Census Web Portal. (https://www.scotlandscensus.gov.uk/webapi/jsf/tableView/tableView.xhtml)(Accessessed: 26 November 2023) Public Health Scotland (2025). Prescriptions in the Community. NHS Scotland Open (https://www.opendata.nhs.scot/dataset/prescriptions-in-the-community)(Accessessed: 26 November 2023) Scottish Government (2020). Scottish Index of Multiple Deprivation 2020v2 data zone lookup file. Scottish Government Supporting Documents. https://www.gov.scot/collections/scottish-index-of-multiple-deprivation-2020/#supportingdocuments(Accessessed UK Government (2025). NHS Health Boards Scotland. Data.gov.uk.https://www.data.gov.uk/dataset/27d0fe5f-79bb-4116-aec9-a8e565ff756a/nhs-health-boards-scotland(Accessessed: 26 November 2023) Rosenthal, M.B. (2004). Seasonal Affective Disorder. In: Encyclopedia of Women’s Health. Springer, Boston, MA. https://doi.org/10.1007/978-0-306-48113-0_391(Accessessed